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Abstract 


Several current challenges in computational fluid dynamics and aerothermodynamics for hypersonic 
vehicle applications are discussed. Example simulations are presented from code validation and code benchmarking 
efforts to illustrate capabilities and limitations. Opportunities to advance the state-of-art in algorithms, grid 
generation and adaptation, and code validation are identified. Highlights of diverse efforts to address these 
challenges are then discussed. One such effort to re-engineer and synthesize the existing analysis capability in 
LAURA, VULCAN, and FUN3D will provide context for these discussions. The critical (and evolving) role of agile 
software engineering practice in the capability enhancement process is also noted. 


Introduction 

Computational fluid dynamics (CFD) is the 
numerical simulation of flowfields through the 
approximate solution of the governing partial 
differential equations for mass, momentum, and 
energy conservation coupled with the appropriate 
relations for thermodynamic and transport properties. 
Aerothermodynamics is the branch of fluid dynamics 
that focuses on the effects of thermodynamic and 
transport models on aerodynamics and heating. It is 
especially focused on conditions of hypersonic 
velocities where the energy content and exchange 
between kinetic, internal, and chemical modes in the 
flow precludes the otherwise common use of calor- 
ically perfect gas assumptions. Computational aero- 
thermodynamics is therefore defined in exactly the 
same manner as CFD, with the added emphasis that 
high temperature gas effects on pressure, skin 
friction, and heat transfer are included in the 
numerical simulation. The fundamental role of 
computational aerothermodynamics is the simulation 
of aerodynamic forces and heating for external and 
internal high speed flows. Reference 1 presents a 
review of recent applications for access to space and 
planetary missions. 

The approximate solution of the governing 
equations requires that the domain be subdivided into 
many small control volumes. The accuracy of the 
solution depends upon a variety of factors; the most 
critical factors are the size of each control volume, 
the orientation of its boundaries relative to a variety 


of flow features, and the order of accuracy of the 
discretization. At the risk of belaboring the obvious, 
note that hypersonic flows routinely involve extremes 
of pressure, density, and temperature separated by 
shocks, expansions, shear layers, and boundary 
layers. These extremes in conditions and topology of 
flow structure complicate computational 
aero thermodynamic simulation. The ability to orient 
a grid with evolving flow structures is particularly 
challenging. For example, accurate simulation of 
heat transfer requires adequate resolution of the 
boundary layer (or merged layer) and accurate 
representation of conditions at the boundary-layer 
edge. Conditions at the boundary-layer edge in turn, 
particularly in the stagnation region of hypersonic 
flows, are dependent on entropy carried along 
streamlines from the shock. Any irregularities in the 
captured shock create associated irregularities in 
entropy that feed the rest of the domain. Further 
complicating the simulation is the fact that the natural 
stability properties of upwind schemes so often used 
in hypersonic applications are in turn most poorly 
realized in the broad, stagnation region of blunt 
bodies where system eigenvalues are small. 

Nevertheless, computational aerothermo- 
dynamic tools that generally do a good job in 
handling the difficulties of steady, laminar, 
hypersonic, blunt body and attached flow simulation 
have evolved over the past 15 years. The greatest 
uncertainty in such applications is the validation of 
physical models for a host of non-equilibrium 
processes 2 . Greater challenges are encountered in the 



simulation of separated flows, complex shock-shock 
and shock - boundary-layer interactions, prediction 
of transition, time-dependent flows, and plasma 
flows. The application of unstructured grid methods 
for such three-dimensional simulations (even the 
relatively simple blunt body problem) has proven to 
be especially challenging when heat transfer is 
required as will be discussed herein. 

The goal of this paper is to reflect on the 
needs of computational aerothermodynamics in order 
to meet these challenges. Some of the approaches 
discussed here are under development within the 
High Energy Flow Solver Synthesis 3 (HEFSS) 
project at NASA Langley Research Center. This 
project combines the physical models and simulation 
capabilities in the structured codes LAURA 4 (focus 
on external, hypersonic flow simulations) and 
VULCAN 5 (focus on internal, scramjet flow 
simulations) into the unstructured flow solver 
FUN3D 6 (perfect gas flow simulations with adjoint- 
based design and optimization capabilities). The 
ultimate goal of this project is to develop robust 
schemes with quantified uncertainties requiring 
minimal user intervention. In this endeavor to re- 
factor existing sources we modularize the code to 
simplify future coupling with other physics. 

Grid Requirements 

Criteria: The HEFSS team is investing in 
unstructured grid technology, including the use of 
mixed elements - tetrahedron (tets), prisms, 
pyramids, and hexahedron (hexes) - as the underlying 
support for future hypersonic flow simulation. We 
defer any investment in patched or overset multiblock 
structured grid systems at this time. The rationale for 
this approach is based on our assessment of current 
and anticipated future capabilities for grid generation 
and grid adaptation. This assessment considers: 

(1) efficiency in the initial generation of a 
“reasonable” grid over complex configurations; 

(2) ability to adapt the grid to flow structures (steady 
and time-dependent) and changes in parameterized 
structures on the vehicles; and (3) solution quality 
and efficiency on a “good” grid using best available 
algorithms. A “reasonable” grid is one that will 
accommodate a converged solution with at least 
approximate resolution of flow structures. It provides 
an initial condition for future adaptation and 
refinement. A “good" grid is one that will produce a 
grid converged solution such that an additional 
refinement will result in dependent variable changes 
within user specified tolerances. 

Initial Generation: Unstructured grid 
generation on complex configurations is already 
judged superior to structured grid capabilities for 
initial grid generation . As an example, consider the 


grid shown in Fig. 1 and Fig. 2 which is discussed in 
Reference 7. A grid around mated ascent vehicles, 
including a fine grid near all solid surfaces to resolve 
boundary layers, was generated in approximately four 
hours (including all user interactions) from a clean 
CAD file. A subjective judgment of the time required 
to create an equivalent, multiblock structured grid 
using continuous or patched interfaces is on the order 
of days to weeks. The time to generate an overset, 
structured grid is estimated to be approximately one 
day. The extra time in these cases derives from the 
topological decisions required by the user so that the 
grid will be sufficiently adaptable to accommodate 
solution quality. 



Fig. 1: GridEx application of unstructured grid for 
simulation of viscous, hypersonic flow over mated 
reusable launch vehicle system. 

Ability to Adapt: The ability to adapt 
structured grids versus unstructured grids for 
hypersonic applications is considered to be a draw at 
present, with unstructured capability advancing 
rapidly. The external hypersonic flow, structured grid 
solver LAURA 4 , for example, routinely applies 
quasi-one-dimensional adaptation as part of the 
solution procedure to align a structured grid with the 
captured bow shock and control near- wall grid size 
and grid stretching factors to obtain grid-converged 
heating rates. However, a lot of time goes into the 




grid generation process in order to enable such 
adaptive capability for complex configurations, and 
some topologies (e.g. elevon gaps as seen in inset of 
Fig. 2) are not amenable to such simple adaptation 
strategies. Many other solution-adaptive tools are 
available to adjust grid spacing based on user 
specified, weighted functions of gradients (feature 
based adaptation) but defining the weight functions is 
often more art than science in order that all flow 
structures are adequately resolved. 


also enforces a surface boundary orthogonality 
condition to improve solution quality of skin friction. 
This issue is likely important for the simulation of 
heating in hypersonic applications as well. 


Mach number - M = 3.0 

m 



Pressure-based adaptation Final grid 
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Fig. 2: Details of unstructured grid in elevon gap for 
RLV application. Scalloped edges are artifacts of 
presentation of tetrahedral facets. 

Feature based adaptation is routinely applied 
to unstructured grids as well. But unstructured grids 
provide additional degrees of freedom that enable 
local enrichment, coarsening, and re-orientation 
(edge-swapping) without propagating these effects 
throughout the domain. However, strategies for 
anisotropic adaptation are required across shocks and 
shear layers to prevent excessive numbers of nodes 
from being drawn into the high gradient regions. As 
an example, the method of Ref. 8 has been applied to 
Mach 3 flow (in viscid) through a scramjet inlet and 
Mach 2 flow (viscous, Re 10000) over a NACA 0012 
airfoil showing very crisp looking captured shocks 
across tightly clustered, high aspect ratio triangles 
(two-dimensional test cases). The algorithm 
examines the eigenvalues and eigenvectors associated 
with the Hessian of each dependent variable and 
derives a transformation based on this geometric 
information that guides directional refinement. It 
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Fig. 3: Mach number contours and feature-based, 
isotropic adapted grid from Venditti 9 for Mach 3, 
inviscid flow over airfoils. 

Adjoint-based grid adaptation is a recent 
development that refines mesh as a function of user 
defined error tolerances. The adjoint equations arise 
from the linearization of the effects of residual 
changes at every node on an output function (lift, 
drag, heating). A user can specify the desired 






tolerance on the output function and the adaptation 
process continues until the tolerance is met or system 
resource limits are exceeded. In this sense, adjoint 
adaptation is an objective error minimization 
technique whereas feature adaptation is a subjective 
error reduction technique. 

An example of adjoint-based adaptation and 
feature-based adaptation applied to supersonic flow 
over a pair of airfoils is shown in Fig. 3 and Fig. 4 
taken from Ref. 9 and 10. The objective function in 
this example is the drag on the lower airfoil. The 
adjoint adaptation yields the same drag coefficient as 
the Hessian based adaptation to about 0.013% using 
nearly a factor of ten fewer nodes. Note that aniso- 
tropic adaptation has not been included in this 
example but the capability is currently being tested 
within HEFSS. 


Adjaint-based adaptation Final grid 

C d = Q.Q766 (lower airjhil) 3810 Nodes 



Fig. 4: Isotropic, adjoint-based grid adaptation to 
obtain grid converged drag on lower airfoil from 
Venditti 9 . Only parts of flow required to get grid 
converged drag on lower airfoil are resolved. 

It is worth noting that the adjoint equation 
set provides capability for configuration design and 
optimization 11 . This capability is currently applied to 
perfect-gas, low speed to transonic applications in 
HEFSS. It has not been applied to hypersonic appli- 
cations because of uncertainty in the choice of 
optimum algorithm (reconstruction, limiting) re- 
quired for flow simulations including strong shocks. 
Still, it is anticipated that this capability will be of 
fundamental value and was one of the major reasons 
for choosing the FUN3D code as the baseline for 
extensions to the hypersonic domain. 


Solution quality and efficiency: 

Structured grid solvers are currently superior to un- 
structured grid solvers as regards quality and 
efficiency for viscous , hypersonic flow applications 
including heat transfer . This assertion, based on 
personal observations, ignores the time required to 
generate an initial grid and assumes that the number 
of nodes in both approaches is comparable. Un- 
structured grid solvers carry an overhead associated 
with connectivity of nodes, edges, and volumes - 
particularly in the case of mixed elements. The con- 
vergence rate of unstructured solvers, especially in 
the presence of strong shocks is slower than 
structured grid algorithms. This behavior is thought 
to be a consequence of greater problems with 
“ringing” associated with the application of limiters 
in the flux reconstruction algorithms and with less 
efficient implicit relaxation algorithms but more 
research is required to confirm such speculation. 

The solution quality concern stems from 
observations of blunt body stagnation region heating 
in 3D using HEFSS on simple test problems 
involving cylinders and spheres. (This issue will be 
discussed in more detail in the Algorithms section.) 
Heating contours in these tests display fair to poor 
symmetry in the stagnation region, depending on the 
reconstruction algorithm, whereas symmetries are 
very good away from the stagnation point. It is 
believed that unstructured-grid-induced irregularities 
in the captured shock promote entropy gradients that 
feed directly to the boundary layer edge in the 
stagnation region. These irregularities tend to be 
dissipated downstream as the entropy layer enters the 
boundary layer. The tested grid in the cylinder case 
was entirely composed of tets derived from an 
adapted, structured grid. The diagonals inserted into 
the structured grid to create the tets were all oriented 
identically; creating a highly biased grid that provides 
an extreme test case. (A literature search for similar 
tests focused on 3D heating was unsuccessful.) 

Caveats: 

(1) The problem of relative motion of two 
rigid bodies where boundary layer resolution is 
required on each body is probably better suited to 
overset grid algorithms. Both structured 12 and 
unstructured 13 grids may be used as the baseline 
discretization. 

(2) The “unstructured” Cartesian grid 
method CART3D 14 is a highly efficient alternative to 
traditional unstructured grid approaches in terms of 
time to generate and refine a grid and in terms of time 
to generate a solution on the grid. This approach 
employs a structured, rectangular grid system 
throughout the domain (Fig. 5) but enables local 
refinements to better resolve regions of high 
curvature in the body and high gradient in the flow. 



Fig. 5: CART3D grid over nose of X-34. 

Preliminary tests of the code for inviscid, 
hypersonic applications indicate it is very attractive 
for initial rapid screening of inviscid aerodynamics 
and can be used to drive boundary-layer codes to 
obtain convective heating using the axisymmetric 
analog, notwithstanding some concerns associated 
with effects of stair-stepping cuts across curved 
surfaces and shocks. The inviscid plus boundary 
layer analysis (structured or unstructured) is still the 
most efficient method for defining heating levels over 
large acreage areas of attached flows over hypersonic 
entry vehicles. 

(3) The perspective throughout this 
presentation is on Reynolds-Averaged Navier-Stokes 
equations which require fine resolution in only one 
direction (anisotropic adaptation to high aspect-ratio 
cells) to resolve boundary-layers, shear layers, and 
shocks (see caveat (4)). External hypersonic flows 
generally feature broad, thin shock layers that also 
are most efficiently resolved with high aspect ratio 
cells (though not as high as required for the 
boundary-layer). Future applications using Large- 
Eddy Simulations 15 (LES) require a more isotropic 
resolution of sub-boundary-layer scale fluctuations 
across the boundary layer for turbulent flow 
solutions. This path is not being pursued now in 
HEFSS because availability of resources for even 
modestly complex configurations is too far (-15 
years) in the future 16 . However, LES for hypersonic 
applications should be investigated on unit, bench- 
mark problems to prepare for that future. 

(4) It is not clear if hexagonal or prismatic 
elements with faces parallel and orthogonal to shocks 
are preferred for shock capturing in the context of an 
otherwise unstructured system. Structured-grid 
algorithm experience with LAURA shows that fine 
grid resolution across the shock is not as important as 


alignment with the shock. Skewed faces across a 
captured shock may introduce kinks with effects that 
are evident downstream. Anisotropic adaptation as 
cited earlier in Ref. 8 appears be an easier approach 
to the challenge at the expense of additional grid 
resources. More focused research on 3D hypersonic, 
blunt body flows with attention to heat transfer is 
required to resolve this issue. 

Algorithm Concerns 

Reconstruction and Limiters: Recon- 
struction refers to the algorithm used to derive the 
right and left states across a cell wall required in the 
flux definition in finite- volume, approximate 
Riemann solvers. Limiters refer to the algorithm used 
to constrain the reconstructed right and left states to 
prevent introduction of new extrema in the flowfield. 
Engagement of a limiter typically reduces the local 
order of accuracy of the reconstructed flux from 2 nd 
to 1 st order. The effects of various linear recon- 
struction techniques on unstructured meshes are 
discussed in Ref. 17 by Aftosmis et al. They note an 
increased degradation of solution quality associated 
with engagement of limiters on triangular (unstruc- 
tured) meshes as compared to quadrilateral (struc- 
tured) meshes. They also note that the use of 
directional limiting (local engagement of the limiter 
for each edge direction) as opposed to a global limit- 
ing (engagement of the limiter considering all 
directions emanating from a node) produced higher 
quality solutions. These observations are consistent 
with behavior observed in the HEFSS benchmark 
simulations, to be discussed subsequently. 

A high level overview of flux definition in 
HEFSS is provided below though several important 
details (eigenvalue limiting, flux reconstruction 
variable set, eigenvector scaling) are omitted. The 
detail here is sufficient to provide context for the 
sensitivity of stagnation heating as a function of the 
reconstruction algorithm. The basic algorithm uses 
the divergence theorem (Green-Gauss) for viscous 
gradients and Roe’s method 18 for inviscid flux: 

f = \ +f„ - R-'|A|(dq - dq_)] (1) 

where 

dq = R(q L - q R ) 

R and R' 1 are the row and column eigenvector 
matrices of the flux Jacobian and A is the cor- 
responding diagonal matrix of eigenvalues. Two 
options are offered here to reconstruct the left (L) and 



right (R) states. The first (Option 1) uses a least 
squares gradient 19 with a limiter function by Barth 20 

or Venkatakrishnan 21 . The term dq lim in Eq. 1 is 
identically zero for Option 1 . The second (Option 2) 
uses the nodal values for the left and right states. 
Option 2 achieves second-order spatial accuracy with 
the Symmetric Total Variation Diminishing 22 
(STVD) scheme applied as a directional limiter 
(function of edge orientation) as follows: 

dqlim = 2 rninmod(dq L , dq R , dq, ^<dq L +dq R )) 
dq M = R[ds • Vq M ] 

ds = [(x L - x R )i + (y L - y R )j + (z L - z R )k] 



The gradients Vq M at nodes (M = L or R) 
are computed using the identical, unlimited least 
squares formulation from Option 1 . The STVD 
limiting in Option 2 and the Venkatakrishnan limiting 
in Option 1 are observed to diminish but not 
eliminate overshoots and undershoots in temperature 
and pressure in the vicinity of strong shocks. A 
pressure ratio weighting function is used to smoothly 
switch to first-order in the vicinity of strong shocks to 
augment the limiting process in Option 2. Options 1 
and 2 use identical formulation of the viscous terms. 



Fig. 6: Structured grid and equivalent uniformly 
biased unstructured grid in symmetry plane over 
cylinder for hypersonic test case. 


Fig. 7: Uniformly biased, unstructured surface grid 
derived from structured grid with ten span wise cells. 

Benchmark Tests: Benchmark cases for 
hypersonic flow over a sphere and a cylinder at 
V„ = 5000 m/s, p„ = 0.001 kg/m 3 , T„ = 200 K 
(approximately Mach 17) consistent with earlier 
FAURA benchmarks, have been intensively worked 
during the HEFSS code development process. Cases 
presented here are for single species N 2 but 
thermochemical equilibrium and nonequilibrium 
options are operational. The structured FAURA grid, 
adapted to align with the bow shock, is shown in the 
bottom part of Fig. 6. The unstructured grid for the 
HEFSS tests (top of Fig. 6) is formed directly from 
the structured grid. Both grids are three-dimensional 
(structured hexahedra (hexes), unstructured 
tetrahedral (tets)). The unstructured surface grid is 
shown in Fig. 7. The grid is highly biased with 
diagonals running in the same direction. The biasing 
is intentionally retained to test simulation quality 
under adverse conditions. Span wise symmetry of the 
solution across the bias is easily monitored. 
Comparisons to the structured grid FAURA results 
are also available in the evaluation process. 

The single spanwise cell simulation results 
with Option 2 (STVD) are shown in Fig. 8. In this 
case there are no interior nodes in the spanwise 
direction. The HEFSS results for heating and 
pressure at surface nodes (indicated by symbol 
location) on the front and back plane are presented; 
the pressure symbols on the front and back plane 
exactly match, the heating levels exhibit a very small 
offset (less than half a symbol height). Good 
agreement with the FAURA reference is achieved 
and the near perfect over plotting of the front and 
back plane results indicates excellent spanwise 
constancy. 



Fig. 8: Pressure and heating distribution for cylinder 
test case with 3D, uniformly biased, unstructured grid 
across one span wise cell. Front and back symmetry 
planes results essentially overplot as expected for 
perfect spanwise uniformity. 
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Fig. 9: Pressure and heating distribution for cylinder 
test case with 3D, uniformly biased, unstructured grid 
across ten spanwise cells with Option 2. Every 
surface node is plotted. 

The ten spanwise cell simulation results with 
Option 2 (STVD) are shown in Fig. 9. In this case 
there are nine interior nodes on each spanwise line as 
shown in Fig. 7. The value at every surface node 
from the HEFSS simulation is plotted. The eleven 
spanwise nodes for pressure again show excellent 
constancy but the heating shows significant spanwise 
variation approaching the y=0 plane (Fig. 10), with 
maximum spread in the stagnation region and 
solution quality improving away from the stagnation 
region. The average value of both the pressure and 
heating continue to be in good agreement with the 
LAURA results. The equivalent result using Option 1 
with the Venkatakrishnan limiter is shown in Fig. 11. 


In this case, pressure constancy is slightly degraded 
and more significant asymmetry in heat transfer is 
evident. Though not shown, the Option 1 results with 
the Barth limiter are of poorer quality and appear 
more like a first-order result. In all cases, the 
symmetry of shock layer contour plots of pressure, 
temperature, and Mach number show consistent 
trends; the spanwise symmetry properties of the 
Option 2 path are superior to all other tested options. 



Fig. 10: Surface heating contours on cylinder from 
Option 2. 



Fig. 11: Same as Fig. 9 with Option 1 and Venkata- 
krishnan limiter. 

These results demonstrate that solution 
quality for stagnation point heating behind strong 
shocks using tets all the way to the wall is poor using 
any of the reconstruction options discussed herein. 
Improvement is expected if prisms are applied in the 
near wall region but this option was not available in 
time to be discussed here within HEFSS. It is unfair 
to lay all of the blame on the Green-Gauss formula- 
tion of viscous gradients across high aspect ratio cells 
when the results are so obviously sensitive to the 


higher order inviscid formulation. In any case, the 
freedom to utilize tets across the entire domain - 
even to the wall - provides the greatest flexibility for 
future grid adaptation. Note that the first-order, 
edge-based reconstruction in Eq. 1 only involves 
information at nodes R and L. A conjecture is offered 
that a truly multi-dimensional reconstruction of the 
flux at every face could overcome some of the 
simulation quality issues noted above. A multi- 
dimensional upwind (MDU) algorithm for high speed 
flows has recently been developed 23,24 ; it is not yet 
known how it would perform on the ten-spanwise- 
cell cylinder test problem for hypersonic, blunt-body 
heating. 

Shock capturing: Special algorithms for 
shock capturing are obviated if anisotropic adaptation 
and/or a multi-dimensional reconstruction algorithm 
provide acceptable solution quality. If locally 
orthogonal grid alignment is required as well (as 
discussed in Caveat (4) above) then discrete shock 
detection and tracking 25 should be considered. 
Ultimately, the algorithmic complexity and overhead 
of adapting suitable grids around shocks may be no 
less onerous than that required to detect and track 
discrete shocks as they move across the background 
grid and modify the reconstruction algorithm when a 
discrete shock cuts an edge. 

Code Validation 

Validation of computational aerothermo- 
dynamics codes requires access to both ground- and 
flight-data. In some cases, Direct Numerical 
Simulation (DNS) of simple, “unit” problems is 
recommended to calibrate physical models. 

Flight data: Comparison to flight data is the 
most important source of validation for the suite of 
physical models required for aerothermodynamic 
simulation. Flight data 1,2 (STS - accelerations and 
heating rates, FIRE 2 - convective and radiative 
heating, Mars Pathfinder - accelerations, RAM-C - 
electron number densities) are important, particularly 
for new vehicles that are modest perturbations of 
tested configurations. However, boundary conditions 
are not always well defined in flight and integrated, 
multi-physics effects usually compromise validation 
of component physical models. For example, in the 
FIRE 2 data set there are multi-temperature, chemical 
kinetics models; thermal relaxation models; and 
surface catalytic boundary conditions (among others) 
that must all be right to agree with measured data. 
Offsetting errors may not be identified - one may 
only conclude that a specific set of models is 
consistent with the measured data set. 


Facility data: Ground-based tests provide 
the best opportunity to measure boundary conditions 
and probe the flowfield with sufficient detail to 
validate physical models. The facilities are not 
perfectly quiet and often the flows are neither 
perfectly uniform nor perfectly steady. As enthalpies 
are increased flow quality and test times become 
more problematic. New facilities may alleviate some 
of these difficulties 26 . However, the problematic 
conditions are accessible to measurement and 
subsequent specification in the simulation. One 
strives for benchmark tests that are simple to simulate 
(axisymmetric, steady, uniform, laminar - if 
turbulence is not specifically being studied) but in 
reality the validation process must accommodate the 
complexities of the test. Consider the following 
validation tests as an example of this process. 

A series of experiments to measure pressure 
and heating for code validation involving hypersonic, 
laminar, separated flows was conducted at the 
Calspan-University at Buffalo Research Center 
(CUBRC) in the Farge Energy National Shock 
(FENS) tunnel. These tests were the focus of a 
“workshop” for hypersonic applications organized by 
Holden and Harvey 27 that should be used as a model 
for ongoing code validation work. (Its format evolved 
from earlier, groundbreaking workshops with broader 
focus on a variety of high-speed simulations 28 .) The 
experimental data were not available to the authors 
until the conclusion of the workshop session, held at 
the AIAA Aerospace Sciences Meeting in Reno 
(2001). Participants transferred their results in pre- 
scribed formats to the organizers for public display in 
the final presentation of the session with respect to 
experimental data and others’ simulations. One set of 
experiments involve Mach 9.5 and Mach 11.3 N 2 
flow over a hollow cylinder-flare with 30° flare angle 
at several Reynolds numbers sustaining laminar, 
separated flow. Truncated and extended flare 
configurations were considered. A second set of 
experiments, at similar conditions, involved flow 
over a sharp double cone (Fig. 12) with fore-cone 
angle of 25° and aft-cone angle of 55°. Both sets of 
experiments included 30° compressions. 

Focation of separation is extremely sensitive 
to the level of grid refinement in the numerical 
predictions. For example, at conditions for 
simulation of Run 28 for the sharp, double cone 
(V,„ = 2664. m/s, p x , = 0.000655 kg/m 3 , T„ = 185.6 K, 
M, x = 9.59, and Re,* = 144010. m’ 1 ), a grid converged 
solution for onset of separation is not achieved until a 
grid density equal to 512 streamwise by 256 normal. 
(A plot of this grid in half column format would 
show almost no white space between grid lines.) A 
simulation using half the number of nodes in both 
coordinate directions yields a 20% smaller separation 



zone. The reflected wave pattern at reattachment 
(Fig. 13) is very complex. The value of an adjoint 
based grid adaptation capability (recall Fig. 4 relative 
to Fig. 3) for a problem like this cannot be overstated. 
The current grid is highly refined even in areas of 
smooth flow and without adaptation, including 
automatic coarsening, the requirements for this same 
problem at angle of attack (3D flow) is estimated to 
be of order of 10 6 nodes. 



Fig. 12: Sharp, double cone (SDC) and computational 
domain bounding pressure contour solution. 
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Fig. 13: Detail of interaction region for Run 28 with 
pressure contours (thin, multi-colored lines), 
streamlines (thin, black lines), and sonic line (thick, 
red line). 
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Fig. 14: Surface pressure and heating distributions for 
Run 28. 

The most interesting development from this 
study concerned the simulation of heat transfer and 
pressure along the forward portion of the sharp cone. 
This region in front of the separation region was 
expected to be the easiest to “get right” and yet was 
consistently overpredicted by the participants (Fig. 
14). A follow-on computational study by Candler et 
al. 29 showed that thermal nonequilibrium effects in 
the expansion from the shock tube nozzle and 
turbulent transition on the tunnel walls were required 
in the simulation to recover proper freestream 
conditions. When these effects were included in the 
simulation with a model for vibrational energy 
accommodation at the surface, the measured 
forebody heating and pressures were recovered. 

These phenomena had less severe impact on the flow 
in the interaction region where simulations had been 
in relatively good agreement with experimental data 
- the thermal effects had sufficient time to equilibrate 
in the recirculating flow. 

These follow-on efforts by Candler et al. 29 
and in a related NATO Research and Technology 
Organization (RTO) study 30 underscore the 
proposition that future code validation tests will 
require more complete simulation of the facility flow 
in addition to simulation around models. Consider for 
example recent comparisons with high enthalpy tests 
in T5 31 and X2 32 . The quality of the comparisons is 
encouraging but uncertainty exists in the 
specification of boundary conditions. Benchmark 
simulations of facility flow (time accurate, three- 
dimensional, appropriate thermochemical models) 
over a range of operating conditions and atmospheres 


will need to be part and parcel of future validation 
studies. As in the case of flight data, there are multi- 
physics effects that obscure the validation; however, 
a richer opportunity to make measurements can be 
exploited to better quantify uncertainties. 

DNS: On top of any list for code validation 
needs is prediction of hypersonic transition and 
turbulent heating levels. This is especially true for 
vehicles that are designed to cruise at hypersonic 
conditions (NASP, X-43) as opposed to enduring a 
hypersonic environment during ascent and entry 
(STS). Understanding the effects of turbulence on 
mixing and ignition delay are critical priorities for the 
scramjet community. The goal here is to calibrate 
the sub-grid scale physics in an LES or multi- 
equation turbulence models using DNS and then 
validate the new model with experiment. In theory, 
DNS can gather statistics on a fine scale inaccessible 
to experimental methods. Dong and Zhong in this 
regard discuss recent efforts 33 for structured grid, 
high-order (5 th ), shock fitted algorithm to simulate 
hypersonic, blunt body flow. It is interesting to note 
that the problems associated with high-order shock 
capturing are avoided in their shock fitted approach. 

Software Engineering 

As computational aerothermodynamic 
simulation codes grow in complexity, issues of 
software engineering (how one develops, integrates, 
and maintains software) become critically important. 
The HEFSS team has struggled with software 
engineering issues that pit code modularity and 
maintainability versus code efficiency. At one 
extreme, future research and development users of 
the software need to be able to modify algorithms, 
add functionality, couple new physics, and expand 
boundary condition options with confidence that the 
task will not take weeks to figure out and their 
changes will not break something else. At the other 
extreme, application users of the code would like for 
it to run as quickly as possible so that high fidelity 
information can be brought more effectively into the 
design process. Our team is far from unique in facing 
these issues. We have adopted some strategies that 
work well for us given our initial conditions (diverse 
expertise in three codes) and environment 
(government lab). 

The HEFSS team includes a dozen people 
across four branches at Langley who actively create 
and modify code on an ongoing basis. It is adopting 
Extreme Programming (XP) methods in this 
endeavor within a research environment 3435 . This 
approach has been a major culture shift - but it is one 
we believe essential in order to engage the many 
elements of this project. For example, the merging of 


mixed-elements, adjoint-based optimization 
algorithms, implicit line solves, and hypersonic 
physical models progressed simultaneously on a 
continuously upgraded suite of modules. Two 
elements of XP that have been critical in achieving 
our goals are continuous testing/builds and pair 
programming. 

Testing: Source code is maintained in a 
Concurrent Version System (CVS - 
http://www.cvshome.org ) repository, (one of several 
revision control systems available for software 
management). Every “expert” on the team is 
responsible for defining a sufficient number of tests 
to cover software functionality within their domain. 
The tests (over 250 and counting) run on a 
continuous basis and check compilation and 
execution for a variety of applications on several 
platforms using a fresh checkout of code from the 
repository. All team members are encouraged to 
“commit” their work to the repository as it develops 
as quickly as possible without any requirement to 
obtain approval from a software manager. (A CVS 
“commit” revises software in the repository accessed 
by all users.) In this environment, one finds that if 
too much time passes between “commits” of new 
material then the resolution of software conflicts 
becomes more tedious and difficult. If “commits” are 
too careless, without adequate checks using the test 
suite, the continuous builds and tests will catch a 
failure and automatically notify the entire team. In 
this regard, peer pressure works to encourage good 
programming practice. 

Pair Programming: Pair programming is 
an effective way of growing knowledge about code. 
Team members work in pairs to develop code and 
implement algorithms in real time. More people 
have detailed knowledge of code and programming 
strategies because of these pair programming 
sessions. The non-expert participant (for a particular 
algorithm or physical model) generally is better 
suited to define the level of comments required to 
instruct future developers. 

Implementation: Our implementation of 
XP arose from a series of seminars organized by two 
advocates of XP within our team. We have developed 
a strategy that works for us and continues to evolve. 
In my opinion, the process would never have been 
effectively implemented had it been mandated from 
above - a healthy group dynamic is critical to its 
success. 

Concluding Remarks 

Opportunities to advance the state-of-art in 
algorithms, grid generation and adaptation, and code 
validation are identified with emphasis on the needs 


of computational aerothermodynamics. Unstructured 
grid algorithms with adjoint based adaptation 
capabilities provide the greatest opportunity for next 
generation hypersonic flow simulation and vehicle 
design. The primary justifications for investment in 
such technology arise from the efficiency in genera- 
ting an initial, workable grid from CAD files (as 
compared to multiblock structured, CO continuous or 
overset grids) and in the ability to adapt these grids 
using adjoint algorithms. Adjoint adaptation provides 
an objective metric of grid quality based on user 
specified tolerances on errors associated with grid 
convergence. It provides logic for coarsening as well 
as refining grids so that only the grid resources 
required for attaining the target error norms are used 
in the simulation. However, unstructured grid 
algorithms in hypersonic applications exhibit short- 
comings in regard to simulation of heating - a 
fundamental requirement of computational aero- 
thermodynamics. These problems stem from 
difficulties in providing high-order accurate 
representation of gradients using high aspect ratio tets 
(usually overcome with use of prisms in the near wall 
region), entropy gradients associated with jagged 
shock captures flowing into the boundary layer, and 
the adverse effects of limiters in an unstructured 
system. 

Examples of algorithmic effects on stag- 
nation region heating are provided herein. Tests are 
conducted on hypersonic flow over cylinders using a 
highly biased, unstructured grid discretization. The 
stagnation region heating is especially sensitive to 
details of the inviscid flux definition. A Symmetric 
Total Variation Diminishing (STVD) algorithm 
provides the best results of options tested in this 
environment (biased tets all the way to the wall); 
however, solution quality (as judged by span wise 
constancy) is marginally acceptable for this case. 
Further improvement using multi-dimensional, 
upwind algorithms and more flexible, mixed-element 
grid topologies will be sought in future work. 
Anisotropic, mixed-element grid adaptation that 
better aligns with shock and high gradient regions 
may relax the need for additional algorithmic 
improvement in these hypersonic applications. 

Validation of new and existing algorithms 
continues to be a critical requirement of computa- 
tional aerothermodynamics. Code validation based on 
experimental data of “unit” flow problems is an ideal 
that enables evaluation of specific simulation 
capabilities and requirements. As facilities strive to 
push the gas harder to achieve high enthalpy 
hypersonic flow, it becomes increasingly difficult to 
provide well defined boundary and initial conditions 
for “simple”, axisymmetric flow simulation. The 
CFD validation process must address the need to 


extend the simulation across the facility itself to 
define the core flow feeding the “simple” configura- 
tion and test. An example of this process by Candler 
et al. 29 pertaining to validation data from Holden et 
al. 27 showed that facility simulation was required to 
better quantify thermal nonequilibrium effects in the 
free stream and their consequent effect on heating 
and pressure on the forward portion of a sharp cone. 

Finally, the critical (and evolving) role of 
agile software engineering practice in the develop- 
ment of simulation capability is noted. Introduction 
of extreme programming (XP) into a research- 
oriented team was initially met with skepticism but 
has been welcomed as experience is gained. The 
primary reason for the change of heart is the 
recognition of how quickly diverse capabilities by 
many authors can be merged in a shared suite of 
modules providing real time access to the latest 
advances in physical models, boundary conditions, 
and algorithms. 

Acknowledgements 

Bill Jones of FaRC provided Fig. 1 and Fig. 
2. David Venditti of MIT provided Fig. 3 and Fig. 4. 
Jim Weilmuenster of FaRC provided Fig. 5. 


References 

1 Gnoffo, P. A.; Weilmuenster, K. J.; Hamilton, H. H. 
II; Olynick, D. R.; Venkatapathy, E.: “Computational 
Aerothermodynamic Design Issues for Hypersonic 
Vehicles,” JSR, Vol. 36, No. 1, pp 21-43, Jan. 1999. 
2 Gnoffo, Peter A.: “Planetary-Entry Gas Dynamics,” 
Ann. Rev. Fluid Mech. 1999, 31:459-94. 

3 Thomas, J. F. et al.: “Opportunities for 
Breakthroughs in Farge-Scale Computational 
Simulation and Design,” NASA TM 2002-211747, 
June 2002. 

4 Cheatwood, F. M. and Gnoffo, P. A.: “User’s 
Manual for the Fangley Aerothermodynamic Upwind 
Relaxation Algorithm (FAURA)”, NASA TM 4674, 
April 1996. 

5 White, J. A.; Morrison, J. H.: “A Pseudo-Temporal 
Multi-Grid Relaxation Scheme for Solving the 
Parabolized Navier-Stokes Equations,” AIAA Paper 
2002-0143, Jan 2001. 

6 Anderson, W. K.; Bonhaus, D. F.: “An Implicit 
Upwind Algorithm for Computing Turbulent Flows 
on Unstructured Grids,” Computers and Fluids, Vol. 
23, No. l,pp 1-21, 1994. 

7 Jones, W. T.: “An Open Framework for 
Unstructured Grid Generation,” AIAA Paper 2002- 
3192, June 2002. 



8 Castro Diaz, M. J.; Hecht, F.; Mohammadi, B.: 
“New Progress in Anisotropic Grid Adaptation for 
Inviscid and Viscous Flows Simulations,” INRIA No. 
2671, Oct. 1995. 

9 Venditti, D. A.: “Grid Adaptation for Functional 
Outputs of Compressible Flow Simulations” PhD 
Dissertation at Massachusetts Institute of 
Technology, Dept, of Aero, and Astro., June 2002. 

10 Venditti, D. A.; Darmofal, D. L.: “Grid Adaptation 
for Functional Outputs: Application to Two- 
Dimensional Inviscid Flows,” J. Comput. Phys., 176, 
pp 40-69, 2002. 

11 Anderson, W. K.; Venkatakrishnan, V.: “Aero- 
dynamic Design Optimization on Unstructured Grids 
with a Continuous Adjoint Formulation,” AIAA 97- 
0643, Jan. 1997. 

12 Atwood, C.: "Computation of a Controlled Store 
Separation from a Cavity," AIAA J. Aircraft, Vol. 32, 
No. 4, July-August 1995, pp. 846-852. 

13 Togashi, F., Nakahashi, K., Ito, Y., Shinbo, Y., 
Iwamiya, T., "Flow Simulation of NAL Experimental 
Supersonic Airplane/Booster Separation Using 
Overset Unstructured Grids", Computers & Fluids, 
Vol. 30, pp. 673-688, 2001. 

14 Aftomis, M.: “3D Applications of a Cartesian Grid 
Euler Method,” AIAA 95-0853, Jan. 1995. See also: 
http://people.nas.nasa.gov/~aftosmis/cart3d/cart3Dho 
me.html 

15 Lesieur, M.; Comte, P.: “Large-Eddy Simulations 
of Compressible Turbulent Flows,” in AGARD 
Report 819, Turbulence in Compressible Flows, pp 4- 
1 -4-39, June 1997. 

16 Lohner, R.: “Generation of Unstructured Grids 
Suitable for RANS Calculations,” AIAA-99-0662, 
Jan. 1999. 

17 Aftosmis, M.; Gaitonde, D.; Tavares, T. S.: 
“Behavior of Linear Reconstruction Techniques on 
Unstructured Meshes,” AIAA J., Vol. 33, No. 11, 
Nov. 1995, pp2038-2049. 

18 Roe, P. L.: “Approximate Riemann Solvers, Para- 
meter Vectors, and Difference Schemes,” J. Comput. 
Phys., vol. 43, no. 2, Oct. 1981, pp. 357-372. 

19 Barth, T. J.: “Numerical Aspects of Computing 
High-Reynolds Number Flows on Unstructured 
Meshes,” AIAA Paper 91-0721, 1991. 

20 Barth, T. J.; Jespersen, D. C.: “The Design and 
Application of Upwind Schemes on Unstructured 
Meshes,” AIAA Paper 89-0366, Jan. 1989. 

21 Venkatakrishnan, V.: “Convergence to Steady- 
State Solutions of the Euler Equations on 
Unstructured Grids with Limiters,” J. Comp. Phys, 
118:120-130,1995. 

22 Yee, H. C.: “Construction of Explicit and Implicit 
Symmetric TVD Schemes and Their Applications,” J. 
Comput. Phys., vol. 68, no. 1, Jan. 1987, pp. 151-179 


23 Van der Weide, E. “Compressible Flow Simulation 
On Unstructured Grids Using Multi-Dimensional 
Upwind Schemes,” Ph.D thesis, Aero. Eng., Delft 
University of Technology, 1998. 

24 Houtman, E.M.; Van der Weide E.; Deconinck, H. 
and Bakker, P.O.: “Computational Analysis of Base 
Flow/Jet Plume Interaction,” Proc. 3d European 
Symp. on Aerothermodynamics for Space 
Vehicles, ESA SP-426, pp 605-612, 1998. 

25 Moretti, G. Efficient Euler Solver with Many 
Applications,” AIAA Journal, Vol 26. No. 6 (1988), 
p655-660. 

26 Mansfield, D. K.; Grinstead, P. J.; Howard, P. J.; 
Brown, G. L.; Girgis, I.; Miles, R. B.; Lipinski, R.; 
Pena, G.; Schneider, L.; Howard, R.: “A 1 MW 
Radiatively-Driven Hypersonic Wind Tunnel 
Experiment,” AIAA 2002-3130, June 2002. 

27 Holden, M. and Harvey, J.: “Comparisons Between 
Experimental Measurements over Cone/Cone and 
Cylinder /Flare Configurations and Predictions 
Employing DSMC and Navier-Stokes Solvers,” 
AIAA 2001-1031, January 2001. 

28 Borrelli, S.; Grasso, F.; Marini, M.; Periaux, J. 
(eds.): First Europe-US High Speed Flow Field 
Database Workshop , Published by AIAA, ISBN 1- 
56347-341-0, 1997. 

29 Candler, G. V.; Nompelis, I.; Druguet, M.-C.; 
Holden, M. S.; Wadhams, T. P.; Boyd, I. D.; Wang, 
W-L.: “CFD validation for hypersonic flight - 
Hypersonic double-cone flow simulations,” AIAA 
Paper 2002-0581, Jan. 2002. 

30 Candler, G. V.; Nompelis, I.; “CFD Validation for 
Hypersonic Flight - Real GasFflows,” AIAA Paper 
2002-0434, Jan. 2002. 

31 Rasheed, A.; Fujii, K.; Hornung, H. G.; Hall, J. L.: 
“Experimental investigation of the flow over a 
toroidal aerocapture ballute,” AIAA Paper 2001- 
2460, June 2001 

32 McIntyre, T. J.; Bishop, A. I.; Rubinsztein-Dunlop, 
H.; Gnoffo, P. A.: “Experimental and Numerical 
Studies of Ionizing Flow in a Super-Orbital 
Expansion Tube,” AIAA 2001-0752, Jan. 2001 

33 Dong, H. and Zhong, X.: “Time-accurate Simula- 
tions of Hypersonic Boundary Layer Stability and 
Transition over Blunt Bodies Using Implicit Parallel 
Algorithms,” AIAA Paper 2002-0156, Jan. 2002. 

34 Kleb, W. L., Nielsen, E. J., Gnoffo, P. A., Park, M. 
A.: “Team Software Development in Support of Fast 
Adaptive AeroSpace Tools (FAAST),” AIAA-2003- 
3978, June 2003. 

35 Wood, W. A.; Kleb, W. L.: “Exploring XP for 
Scientific Research,” IEEE Software, Vol. 20, No. 3, 
pp. 30-36, May/June 2003. 



